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We present a general theoretical formulation, based on 
nonequilibrium Green's functions, for nonlinear DC transport 
in multi-probe mesoscopic conductors. The theory is gauge 
invariant and is useful for the predictions of current-voltage 
characteristics and the nonequilibrium charge pile-ups inside 
the conductor. We have provided a detailed comparison be- 
tween the gauge invariant scattering matrix theory and our 
theory. We have also given several examples where the I-V 
curve can be obtained analytically. The effects of exchange 
and correlation have been considered explicitly. 
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I. INTRODUCTION 

Many practical electronical devices, such as diodes and 
transistors, operate on nonlinear current-voltage (I-V) 
characteristics: /q, — /q({V"/3}), where subscripts a and 
(i denote leads which connect the device to the outside 
world. Classically one can predict the I-V curves by solv- 
ing the coupled equations of classical electron motion 
such as the Boltzman equation, and the Poisson equation 
for the electrostatic potential of the conductor, subject- 
ing to the boundary conditions that at the asymptotic 
region of the lead (3, the external bias voltage is fixed 
at Vg. For coherent quantum conductors in the meso- 
scopic regime, one still must solve the coupled equations 
but the electrons are now quantum entities. Clearly the 
prediction of I-V curves becomes much more difficult in 
the quantum situation. As a consequence, most theoret- 
ical analysis of quantum transport in coherent quantum 
devices do not predict I-V curves: they focus on the lin- 
ear DC conductance Gap whichjcan be calculated from a 
variety of theoretical methodsEla. This is then compared 
with experiments which extract Gap from the measured 
I-V curves at a vanishing bias voltage: la = X]/3 GapVp. 

Because linear conductance Gap does not give the 
whole picture concerning a nonlinear device, it is very im- 
portant to, theoretically, go beyond the linear transport 
regime and predict the whole I-V curve. At nonlinear 
situations one must worry about a fundamental physics 
requirement: the gauge invariant condition which dic- 
tates that the predicted electric current should be the 
same when potential everywhere is ^bufted by a constant 
amounta. Biittiker and co-workeraju have developed a 
scattering matrix theory (SMT) which satisfies the gauge 
invariant condition and predicted the second order non- 
linear conductaace. When the scattering matrix takes a 
specially simples form, e.g. the Breit-Wigner form, the 






full nonljjiear I-V characteristics is also obtaineda. The 
key ideaa'Q, in order to maintain gauge invariance, is to 
include the internal potential landscape into the anal; 
sis which is a result of long range Coulomb interaction: 
Recently this SMT has been extenckd to predict higher 
order weakly nonlinear conductancesO and connections to 
the frameworkjjof response theory has bgen formalizedE. 
In both SMTQ'B and the response theor}£l, one calculates 
the nonlinear conductance perturbatively order by order 
in voltage. Such theories make sense for weakly nonlin- 
ear situations where the external bias is finite but small. 
Hence practically it is very difficult to compute TV curves 
for more general situations. 

From the j-. j— . nonequilib- 

rium Green's functions (NEGF)B~0, Ref. ^ provided an 
analysis of I-V curve in the wide-band limit, where gauge 
invariance was satisfied by including the internal poten- 
tial phenomenologically through the use of a capacitive 
charging model. The charging model is, however, not 
fully nonlinear since the internal potential is treated lin- 
early in it. It is thus important and attractive to further 
develop the NEGF to the fully nonlinear regime for the 
purpose of predicting I-V characteristics of multi-probe 
coherent quantum conductors. It is the purpose of this 
work to provide such a development: we have formulated 
a general gauge invariant nonlinear DC theory based on 
NEGF by treating the nonlinear internal potential from 
first principles. Our theory also goes beyond the wide- 
band limit and can directly predicjt-the I-V curves and 
the nonequilibrium charge pile-upsuEj inside the conduc- 
tor. The former can be expanded to obtain weakly non- 
linear conductances which we shall compare with results 
obtainable from SMT at lower orderalU; while the latter 
gives a voltage dependent nonlinear capacitance coeffi- 
cient which is experimentally measurable. Our theory 
provides a solid base for further numerical predictions 
of I-V curves for complicated coherent quantum device 
geometries. We will provide the theoretical formalism in 
section II. In section HI, we will give detailed comparison 
between the result of our theory and that of SMT. We 
will calculate the TV curve for a number of examples. 
The summary will be given in section IV. 



II. GAUGE INVARIANT FORMULATION 

Let's consider a quantum coherent multi-probe con- 
ductor with the Hamiltonian 

H = ^^tka{t)c\^Cka + Hcen{dn,dV\ 



ka,n 
^0 



■a,nCj^^U"n 



c.c] 



(1) 



G'{E,U) = 



where eka{t) = e^ + qVa- The first term of Eq.(n|) de- 
scribes the probes where DC signal is applied far from 
the conductor; the second term is the general Hamil- 
tonian for the scattering region which is a polynomial 
in {dJjjfLJ- that commutes with the electron number 
operatoillj TV = J^n^li^'n', the last term gives the cou- 
pling between probes and the scattering region with the 
coupling matrix Tka.n- Here c|.^ {cka) is the creation 
(annihilation) operator of electrons inside the a-probe. 
Similarly dj^ (cLJjJs the operator for the scattering re- 
gion. As usualHEj, we assume that at time t — > — cx), 
the system is disconnected. 

The electric current can bcrwipi-tten in terms of Green 
functions in the usual mannertS^ta {^ — 1); 

Mt) = -q f dhTr [G'-{t,h) S<(ii,t) 

+G<(i,ii)S°(ii,i)-I]<(i,ti)G°(ii,i) 
-S;(i,ti)G<(ii,i)] . (2) 

Here the capital G's are various Green's functions (see 
below) and the self-energies are defined as follows: 

^Qmn (^1: ^2) = E ^fcam^fcan5feQ'^'''(il> ^2) (3) 

k 

where 

9Litut2)=if{4)eM~^ [ ^kadt) (4) 

Jt2 

and 



fffca (^1- ^2) = TiO{±h T i2) exp(-^ / e^dt). (5) 

Using Fourier transforms of the various quantities on 
the right hand side of Eq. (||), it is straightforward to 
obtain, 



Ja = ~iqj {dE/2^)Tr{T^{E - gK.)[(G''(^, U) 
~G-{E,U))f{E~qV^)+G<{E)]} . (6) 

where G^{E,U) is the retarded Green's function with 
U = U{r) the electro-static potential build-up inside o 
conductor. In the Hartree approximation it is given b; 



G''{E,U) = 



1 



H -qU -^"^ 



(7) 



where S^ = Ea^^l^ " 9^a) and T^{E) = 
— 2/to[E^(_E)] is the liMwidth function. Within the den- 
sity functional theorynj, we can further include the ex- 
change and correlation effect, 



E-H~qU-V^, 



^r 



(8) 



where Vxc — 5Exc/5p is the potential due to the exchange 
and correlation energy Exc and p is the charge density. It 
is worth to emphasize that a most important denaijture 
of our theory from the previous NEGF analysialjO'Ej is 
that we explicitly include the internal potential landscape 
C/(r) into the Green's functions self-consistently. This is 
the crucial step in the development of a gauge invariant 
nonlinear DC theory. At Hartree level U{t) is determined 
by the self-consistent Poisson equation 



V2[/ = -47rp = Amq / [dE / 2tt)G< [E , U) 



(9) 



Within Hartree approximation, G^ is related to the re- 
tarded and advanced Green's functions G*" and G°, 

G<{E,U) = G'Y.'^Tp{E-qVp)f{E-qVp)G'' . (10) 

P 

Eq. (|9|) is, in general, a nonlinear equation because G*"'" 
depends on t/(r) (see Eq.(0)). Using Eq. ( p^ ) we reduce 
the current into the following form, 

Jc^-qY. j{dE/2^)Tr[V^G^TpG'%f^ - fp) (11) 

where we have used the notation Ta = Ta{E — qVa), 
r ~ J2a ^a, and fp = f{E — qVp). To make connection 
with the SMT, we introduce the screened transmission 
function 

^a/3 = Tr[r^G'\Tdo.p - r^)G"]/(4^2) (^3) 

we tben arrive at the familiar form of the current in 
SMlti, 



Jq = -27r9^ / dEfpAap 



(13) 



Eqs. ([TllJTpl) completely determines the nonlinear I- 
V characteristics of an arbitrary multi-probe conductor, 
they form the basic equations of our theory. The self- 
consistent nature of the problem is clear: one must solve 
the quantum scattering problem (the Green's functions) 
in conjunction with the Poisson equation. It is easy to 
prove that the current expression Eq.(|ll|) is gauge invari- 
ant: shifting the potential everywhere by a constant V ^ 
U ^ U + V &iYdVa ^Va + V, Ja from Eq. (|llj) re- 
mains the same. Note in Eq.([ll|) the quantity T depends 
on voltage and without such a voltage dependence, the 
gauge invariance can not be satisfied. On a technical side, 
Eqs. (|ll|,0,^ also form a basis for numerical predictions 
of I-V curves. For instance one can compute the vari- 
ous Green's functions G and the coupling matrix_r for 
multi-probe conductors using tight-binding modelsa; and 
the Poisson equation-ean be solved using very powerful 
numerical technique£2l. 



III. RESULTS 

The main thrust of the previous section (and of this 
work) is the solution of the gauge invariance problem 
for DC nonlinear transport in general terms of bias volt- 
ages. In the rest of the paper we shall derive analyti- 
cal expressions for a number of examples which can be 
solved inxlosed form. At low bias where SMIH and linear 
responseQ are applicable for weakly nonlinear situation, 
we show that our general formula reduces and becomes 
compatible with them. But for higher bias where these 
previous theories are not applicable, our theory becomes 
an unique approach for analyzing nonlinear DC quantum 
transport. Hence we derive the general current-voltage 
characteristics for the entire range of nonlinearity for a 
tunneling device, and prove that results obtained with 
or without gauge invariance can differ substantially not 
only quantitatively but also qualitatively. 



A. Weakly nonlinear regime 

For weak nonlinearity we can series expand all ijuan- 
tities in terms of the small external bias voltageS and 
obtain results order by jOrder: this is precisehf the ap- 
proach adapted in SMTeI and response theoryQ. In this 
subsection we confirm that our nonlinear theory indeed 
reduces to these previous approaches at the weakly non- 
linear regime where they are applicable. In particular we 
shall derive analytical expressions for the local density 
of states (LDOS) and the second order weakly nonlinear 
DC conductance, which are the two interesting quantities 
for weakly nonligear regime. j-. 

In both SMTH and response theorjnD, LDOS plays a 
very important role. From our NEGF theory LDOS can 
be easily derived from the right hand side of Eq. JM} , 
which is the charge density, with the help of Eq. ([lO|). 
Here we shall present the explicit expression at the lowest 
orderQ expansion in the external voltage. Hence we seek 
the solution of U{r) in the following form, 



U ^Ueq+^ UaVa + ■;^ X! ""/3^a^/3 






(14) 



where Ueq is the equihbrium PoteaUaJ and Ua{r), MQ^..(r 
are the characteristic potential»LlEl. 



that 
rulei 



and 



It can be shown 
.characteristic potential satisfy the following sum 



y^Mg = 1 



7e/3 



»{/3}i 



(15) 



(16) 



Here the subscript {j3}i is a short notation of I indices 
^,S,r], ■ ■ ■. Expanding G< of Eq. (fih in power series of 



Vq., we can derive the equations for all the characteristic 
potentials. In particular the expansions are facilitated 
by the Dyson equation to the appropriate order (in the 
absence of the exchange and correlation effect): 



G'' = G^o + Go ilU - qU, 



eq 



^K(£;-gF„)-S^(i?)]jGS 

a / 



(17) 



with Gq the equilibrium retarded Green's function, i.e.. 
when U = Ueq. At the lowest order, we thus obtain 



- V^Uo^ir) = 2g2 / dEf{G'oUa,Gl,TG'^ + c.c.)rr 

- 2g2 j dEfiGlodEKG'o^G^o + c.c-)rr 

^ 2g2 j dE[Gl{T^dEf + a£r„/)Gg]„ (18) 

where Fq, and / no longer depend on voltage after the 
expansion; and c.c. denotes complex conjugate. 

The first term on the right hand side of Eq. 
(|l8|), which depends on internal potential u^, de- 
scribes the induced charge density in the conductor. 
It can be simplified using the fact iT = (GJj)"^ — 
(Gg)^^, hence it becomes ^-Ej£S^^,Wrr'Ua{r') where 
H is the Lindhard functiontju defined as !{„> — 
—i f(d£'/27r)/(Gn^^,Gn^,^-Gn^^,Gn^,^). The second and 
third term of ( |l8| ) which do not depend on characteristic 
potential correspond to the charge density due to exter- 
nal injection. They are the local partial density of states 
(LPDOS) dria (r)/di<^|j:alled injectivity in the language of 
the scattering matrixO, 

dne,{x)/dE = -J{dE/27r)[G'[,r^dEfG^] 



{dE/2Tr)f[G'odETo.G^o 
= - /(dii;/27r)/[GS(GSF„ +F„G[5)GS] 



(d£;/27r)/[GS(c'£E'-GSF„ + F^Ggai^S" 
- a^E^GSF - TG^dEK)G^o] ■ (19) 

Comparing fthis result with that derived by SMTeI, the 
SMT resultQ corresponds to the first term on the right 
hand side of Eq. ( |l9| ) . Hence the local partial density of 
states obtained from our general|-|theory is slightly dif- 
ferent from that defined in SMIB. However it can be 
proven that the difference, e.g. the second integral of 
Eq. (pj|), becomes negligible in the limit of large scatter- 
ing volume. The proof follows the approach detailed in 
our earlier work Ref. E4l 

We obtain LDOS dn{r)/dE from Eq. (|l|): 



dn{v)/dE = ^dna/dE 

'(di?/2^)/[G5(GSr + rGg)GS] 
= i /(d£;/27r)/(GSG5-GSGS). (20) 



This is exactly the same as the LDOS we obtained from 
the response theoryQ. The agreement is actuaUy not sur- 
prising: because we are deahng with DC transport where 
there is time reversal symmetry, the result of NEGF 
should be the same as that of the linear response. Ouc-ce.- 
sult of LDOS exactly satisfies the general relationshipOB: 

J2,, nr,r' = dn{v)/dE. 

As a second comparison to the results obtainable from 
SMT and response theory, we now derive a formula for 
the second order nonlinear conductance Gap-^ from our 
theory of Section II. In weakly nonlinear regime, the elec- 
tric current can be expanded in a series form in terms of 
external bias voltages, 

Jq — 2_^ GapVp -f 2_^ GafijVpVj + 2_^ GaPjsVpVjVs + ... 

p Pi Pl& 

(21) 

To compute the second order nonlinear conductance 
Gap-f, it is enough to calculate u^ for the internal po- 
tential. Expanding the general expression for the current 
Eq.(|ll|) to the second order in voltage, it is straightfor- 
ward to obtain the second order nonlinear conductance: 

Go.p^ = -q" fidE/2TT)Tr [9i=;r„GS(W„7 - rj)G^dEfSo.p 

- r„GS(M^ - dE^)G'^{TSo,-, - r^)G^odEf 

+ ^aGQds^jGQdEfSap — TaGQdE^jGQdEfSpj 

+ (i/2)r„GS(r<5„^ - r^)G^odlfSp^ 

- r„G5(r5„^ - r^)GS(u0 - ^E^)G^^Ef] (22) 

This result is gauge invariant as one can explicitly, 
verify that it satisfies the gauge invariant conditiona 
J2pi^aP7 + C!a^) — 0. This result agrees with that de- 



IT^i 



rived from SMTq if we neglect the terms involving Oe^ 
and 9_eS. In that case, we obtain 



show that the nonlinear characteristic potentials Uap... 
satisfy a Poisson equation similar to (y) with different 
source terms dnap.../dE which correspond to nonlinear 
LDOS. The higher order nonlinear coefficient Gap... can 
be obtained in similar fashion as we have done here for 
the second order coefficient. 



B. I-V curve in the w^ideband limit at zero 
temperature 

In the last subsection we examined the limit of weak 
nonlinearity. However the main advance we have ob- 
tained from the gauge-invariant NEGF formalism devel- 
oped in Section II is to be able to predict the full non- 
linear current-voltage (I-V) curves. Several analysis will 
be presented in this and the next subsections for the I-V 
curves. j— . 

In the commonly used wideband limitll3, the coupling 
matrix F is independent of energy which drastically sim- 
plifies the algebra. The wideband limit corresponds to 
cases where the probes have no feature, thus the inter- 
nal potential U{r) becomes just a space- independent con- 
stant Uo (the value of Uo depends on the voltages {14}). 
In wideband limit the steady state Green's function takes 
a very simple form, Gq = 1/{E — Eq + JF/2), thus the 
integral in Eq. (|ll|) can be done exactly. We obtain. 



(7q — 



p 



F^)Fq, arctan 



AE - qUo + qVp 



F/2 



(24) 



where AE = Ee — E^ and the constant Uo 
is determined by the charge conservation condition 
fdETr[G<{E,V)] ^ J dETr[G<{E,0)], i.e. 



y Tp arctan 

P 



AE - qUo + qVp 



To 



r/2 



obtain 



= F arctan 



A^ 



F/2 



(25) 



this 



Gap^ = q' j {dE/2^)Tr[{GlTaGlupGl + Ggw^^GgF^GS 
l/2GQFQGoGQ(5/3-y — 1/2GqGoFqGq(5/3.^)(F(5q^ — T^)\dEf 



= 2Trq^ / dE{-dEf)[l/2qdEAapSp^ + dv.A^^] 



(23) 



which agrees exactly with the result in Ref. ^ The 
second order non-linear conductance has been investi- 
gated numerically for several systems using the scatter- 
ing approaches. For detailed discussion of the numerical 
technique needed for the calculation, see ref. |2^. 

Finally, we comment that by expanding Eq.(M to 
higher order in terms of voltage, it is straightforward to 



equation, the quasi-neutrality approximationO is assumed 
which neglects the charge polarization in the system in 
addition to the use of total charge neutrality. The gauge 
invariant condition in Eq.(24) is clearly satisfied: raising 
both Vp and Uo by the same amount does not alter the 
current. Eqs.(E4, 03) have been obtained beforea from 
SMT where a Breit-Wigner form of the scattering ma- 
trix is assumed. Hence we may conclude that in this 
sense the wideband limit in NEGF approach is equiva- 
lent to the Breit-Wigner approximation in the scattering 
matrix approach. It is however different from that de- 
rived in Ref. 15^ where a linear charging model is used 
for the internal potential build up. It is not difficult to 
confirm that the result of Ref. ITsI is recovered if we solve 
for Uo in Eq.(E5h to the first order in voltage V, i.e. we 



compute the internal potential Uo by neglecting the con- 
tributions of higher order characteristic potentials Uap.... 
In this limit we obtain Uo — ^ci(^a/^)Va- Substitute 
this into Eq.(24) we arrive at the result of Ref. |l^. This 
exercise also allows us to identify the phcnomcnological 
parameter Ci of Ref. |l^ to be F^, and it indicates that 
the linear charging model for the internal potential is not 
complete for the full nonlinear I-V curve predictions. 

Next let's derive the full nonlinear I-V curve for a quan- 
tum dot with two resonant levels. For two resonant levels 
in a quantum dot the retarded Green's functionCJ can be 
derived to have the following expression: 

r^r 

" [1/{E -El- qU) + 1/{E -E2- qU)]-^ + ir/2 

(26) 

where Ei and E2 are energies of the two resonant levels. 
From Eq.(O), we obtain the current is 

/m[ln(a^-6^) + |ln(^)] 
^ 2,0 ap — 



(27) 



where ap = Ep + qVp - qUo - {Ei + E2)/2 + iT/2 and 
b^ = {El - i?2)^/4 - F2/4. In the quasi-neutrality ap- 
proximation we derive that the internal potential Uo is 
determined by following equation, 



iT ap + b 



Er.Mln(«|-6^) + ^ln( 

= F/m[ln(a2 - 6^) + !£ ln(^)] 



)] 



2b 



-b' 



(28) 



where a ^ Ep - {Ei+ E2)l2 + iV 12. 

Fig.(l) plots the predicted I-V curve Eq.(|2^) with the 
parameter Fi = r2 = 0.1, Ei = 0.3, E2 = 1.4, and 
Ep = 0.11. To compare the I-V curves with and with- 
out gauge invariance, in Fig.(l) we have plotted three 
curves. The dot-dashed line represents the gauge invari- 
ant solution Eq.(|27|). Since this solution is gauge invari- 
ant, we choose V\ —V and V2 = in Fig.(l). Both the 
other two lines (solid and dotted) are for [/q = thus no 
internal potential is taken into account self-consistently. 
Two observations warrant to be discussed. First, the two 
non-self-consistent I-V curves, solid line with Vi = V and 
F2 = and dotted hne with Vi = V/2 and V2 = -V/2, 
give different I-V curves. This is clearly wrong because 
electric current must only depend on the bias voltage 
difference which is V for both curves, and not on the 
choice of the reference point for potential. This is a di- 
rect consequence of the flaw of a non-self-consistent the- 
ory. Second, the qualitative behavior of the current ver- 
sus voltage curve is different. Both curves of the non- 
self-consistent analysis show quantized steps, which the 
self-consistent analysis with quasi-neutrality approxima- 
tion does not give. This difference in qualitative behavior 



can be understood as due to the quasi-neutrality approx- 
imation. The electric current for the incident electrons 
with energy Ep will increase sharply when the chemical 
potential /i = Ep + qV is close to the first resonant en- 
ergy level El. When the internal potential build-up is 
not included in the non-self-consistent analysis, this cur- 
rent saturates after crossing the first level but increases 
again when fi is near the second resonant level. This is 
actually a reasonable picture. In the self-consistent solu- 
tion within quasi-neutrality approximation, however, the 
internal energy Uo solved from quasi-neutrality condition 
Eq. ( p8| ) increases linearly V = kUo with the voltage with 
coefficient k close to one. Hence as the chemical poten- 
tial rises, the resonant level also increases with approx- 
imately the same amount. Thus the second resonance 
level will not be reached (within the quasi-neutrality ap- 
proximation) for the range of voltages of Fig.(l). Finally, 
In Fig. (2), we have plotted the differential conductance 
dl/dV of two non-self-consistent solutions. For the case 
Vi — V and V2 = (solid line), two peaks show up near 

V = 0.2 and V = 1.3. Since the Fermi level in the equi- 
librium is Ep — 0.11, those two peaks reflect the resonant 
behavior when the chemical potential qVi -\- Ep is in line 
with two resonant levels Ei = 0.3 and E2 = 1.4. When 
Vi = V/2 and V2 = —V/2 (dotted line), the chemical po- 
tential is again qVi + Ep = qV/2 + Ep, so we found two 
peaks at T^ = 0.38 and V — 2.6. However, the spacing 
between two peaks are different for two different choice 
of voltage Vi and V2: it is therefore important to include 
the Coulomb interaction so that the theory is gauge in- 
variant. 

Our results strongly suggest that as far as the I- 

V curve prediction is concerned, without interaction 
one violates gauge invariance, but including interaction 
within the quasi-neutrality approximation is still not 
enough as it misses the expected resonance levels in 
the I-V curve- which are often observed in experimen- 
tal situationgH3. Hence it maybe necessary to go beyond 
the quasi-neutrality approximation. 



C. Exchange and correlation effect 

So far the electron-electron interaction which gives rise 
to the internal potential build-up has been treated within 
the Hartree approximation with the quasi-neutrality con- 
dition. In this subsection we examine the effects of ex- 
change and correlation to the nonlinear I-V curves within 
the wideband limit for a resonant tunneling structure. 
We must also go beyond the quasi-neutrality approxima- 
tion. 

Going beyond quasi-neutrality approximation means 
that we must consider the local charge distribution under 
the condition of overall charge neutrality. For a double- 
barrier tunneling structure, let's introduce capacitance 
coefficients Ci and C2 for the left and the right barrier, 
respectively. Then the charge in the quantum well due 



to Coulomb interaction is given b; 



,30 



AQ = -I f{dE/2Tr)[G<{E, Uo) - G<] 

^Ci{Uo-Vi)+C2{Uo-V2) 



(29) 



where AQ is the total charge in the well, Uo is the overall 
shift of the band bottom of the well due to the Coulomb 
interaction, and Gq is the equilibrium lesser Green's 
function. For the system with only one resonant level, 
this equation reduces to 



y Tfj arctan 

/9 



Ef-Ho + qVp 



r/2 



r arctan 



Ep — Hf) 



r/2 



[Ci{Uo - Vi) + C2{Uo - V2)]^T/q 



(30) 



where the wideband limit is assumed and TJo = Eq 
qUo + qVxc- The current is given by 



J a 



-y 

/3 



{TSap — Tp)Ta arctan 



Ep — Hf) 



qVp 



r/2 



(31) 



To plot the I-V curve determined by Eqs.(pO) and (31), 
we use Vxc —t-^Av^'^^Q^^^ i'^ ll^s Hartree-Fock-Slater 
approximationEJta, where 2/3 < a < 1. Parametrized 
by the coupling constants Ti and the geometrical capaci- 
tance coefficients Ci , we thus can calculate the nonlinear 
I-V curves from Eqs.(pO|) and (pl|). Fig. (3) presents cur- 
rent as a function of voltage difference V = 14 — V2. 
We have compared the cases with or without exchange 
and correlation potential Vxc for two different set of sys- 
tem parameters: symmetric barrier with Fi = r2 = 0.5, 
Ci = G2 = 5.0, A£' = -1.0, a = 0.7 (dotted line with 
Vxc, solid line without); and asymmetric barrier with 
Fi = 0.1, F2 = 0.5, Gi = 1.0, G2 = 5.0, A£; = -1.0, 
a = 0.7 (dashed line with Vxc and dot-dashed line with- 
out). We observe that the current for the symmetric 
barrier is much larger than that of the asymmetric bar- 
rier for y > 2, but it can be smaller for smaller bias be- 
tween 1 < y < 2, and becomes larger again at very small 
V < 1. Without the internal potential build-up taken 
into account, it is well known that symmetrical tunnel- 
ing barriers have larger transmission coefficients hence 
larger current than those of asymmetrical barriers. The 
behavior of the FV curves in Fig. (3) at the T^ ^ limit 
is consistent with this picture. However at larger volt- 
ages this expectation may or may not be true, due to the 
nonlinear effects and the internal potential build-up. We 
found that for both symmetrical and asymmetrical bar- 
riers, the exchange and correlation effects are to increase 
the electric current. This is seen more clearly from the 
differential conductance dl/dV versus voltage in Fig. (4). 
Since the exchange and correlation term Vxc is to lower 
the bottom of the conduction band, the peak of dl/dv 
shifts to the small voltage as a result. 

When exchange and correlation effects are included in 
the two-level tunneling system, the Eqs.(p7|) and (pa) are 



modified in a similar way as Eqs.(31) and (|30|). Fig. (5) 
shows the FV curve for two-level system for Fi = F2 = 
0.1, Gi = G2 = 1.0, El = 0.3, E2 = 1.4, Ef = 0.11, and 
a = 0.7. The FV curves with (dotted line) and with- 
out (solid line) Vxc are plotted for comparison, but both 
FV curves now show two steps reflecting the two reso- 
nance levels for tunneling. Hence by going beyond the 
quasi-neutrality approximation, the gauge invariant the- 
ory developed in Section H predicts a "quantized" I-V 
curve which, as mentioned above, is physically reason- 
able. Again, when Vxc is present, the current increases. 
In Fig. (6), we show the differential conductance dl/dV 
for the same system parameters as that of Fig. (5). When 
the exchange and correlation potential is included (dot- 
ted line), it is surprising to observe that there are three 
peaks in dl /dV instead of two (solid line without Vxc)- 
In addition the peaks of dl/dV are shifted towards 
smaller values of bias as compared to the dl/dV curve 
for Vxc — 0. The entire behavior of dl/dV can be under- 
stood as the following. The internal Coulomb potential 
Uo tends to move the resonant level up (to higher energy) 
and the exchange and correlation potential Vxc tends to 
lower it down, these two effects give compensating fac- 
tors to move the resonance levels inside the quantum well. 
At very small voltage, \Vxc\ increases much faster than 
the internal potential Uo does as the bias is increased 
(see inset of Fig. (6)), thus the resonant level Ei moves 
downwards in energy from its "bare" value Ei — 0.3: 
a resonance peak is expected when it is lowered to the 
chemical potential. As the voltage increases such that 
the level Ei is lowered to below the chemical potential, 
dl/dV decreases from the resonance peak consistent with 
the fact of going off resonance. When the voltage in- 
creases further, the resonant levels Ei and E2 start to 
move upward in energy due to the effect of Uo, and it 
will again pass the chemical potential resulting to the 
resonance peak near V = 0.5. Finally when the voltage 
is around V = 2.0, the chemical potential is near the sec- 
ond resonant level and we see a sharp increase of dl/dV 
and the third peak shows up as a result. We conclude 
that the quasi-neutrality condition may need to be ex- 
tended in predicting I-V curves when the system is in the 
tunneling regime. Here we have used a phenomenologi- 
cal but nonlinear capacitance charging model to include 
the charge polarization effects, which are seen to play 
an important role in predicting quantized I-V curves. It 
is further found that the exchange-correlation potential 
Vxc can be quite important as it provides a compensating 
effect to the Hartree internal potential Uo- 



D. I-V curve with Hubbard U term 

To maintain the gauge invariance, we have so far con- 
sidered the Coulomb interaction in the Hartree approxi- 
mation with or without an exchange-correlation term. In 
this subsection, we will consider the on-site Coulomb in- 



teraction in terms of Hubbard U model with the following 
Hamiltonian for Hcen in Eq-djJ): 



Hcen = ^ Eodldcy + Uin^ni 



(32) 



Here we assume that the quantum dot contains one en- 
ergy level Eq with Coulomb repulsion energy Ui which 
accounts for the interaction between different spins. In 
addition to Ui we assume that the long range Coulomb 
potential U between different sites gives an overall a con- 
stant shift Uo to the bottom of conduction band. This is 
similar in spirit to the energy shift A introduced in Ref. 
[l8| , p9| . However, in our case, Uo has to be determined 
self-consistently. 

The current is still determined by Efl.(p^. But now 
the lesser Green's function is given bytHl 



/3 



(33) 



whereBS 

g;(s,c/) 



< n^ > 



E-Eo-Ui-qU + ir/2 

I- < ria > 
E-Eo-qU + iT/2 



and 



< ria- >- 



z / {dE/2n)G< 



(34) 



(35) 



The internal Coulomb potential U can be determined in 
terms of the geometrical capacitances Gi and G2, 

j{dE/2^) ^ G< {E, U) - jidE/2',) ^ G<o - 

<7 <T 

Ci{Uo ~ Vi) + C2{Uo ~ V2) (36) 

from Eg. (|34|) , this condition becomes 

'AE-Ui -qU + qVp' 



e^^j: 



< ria > arctan 



13 

ry < n^ > arctan 



arctan 



r/2 

AE~qU + qVfj 



r/2 



— r 2^(1— < ng >) arctan 



AE -Ui 

rj2 

AE 



r/2 



[C\{Uo - V^i) + C\{Uo - Vi)]27rr/g 



(37) 



With the potential U determined this way, the current is 
finally written as 



Jq, — 



ttF 



^(^Safl -Tfj)Ta 



(y < rig > arctan 



2^(1— < ng >) arctan 



AE-Ui-qU + qVp 
Tj2 

'AE-Ui+qVi 



r/2 



) (38) 



which can be calculated numerically using Eqs.(p3), (p^), 
and (p8|). In Fig. (7), we have plotted the differential 
conductance versus voltage for Fi = r2 = 0.5, Gi — 
G2 == 5, AE = -1.0, and Ui = 4.0. As expected, there 
are two peaks corresponding to two different energies Eq 
and Eq + Ui. Since the Coulomb interaction U increases 
linearly with voltage 2U w V, the separation between 
the two peaks becomes 2Ui. 



IV. SUMMARY 

In this work we have developed a general gauge invari- 
ant nonlinear DC transport theory based on the nonequi- 
librium Green's functions. The main idea of this devel- 
opment is to self-consistently couple the NEGF with the 
proper Poisson equation for the internal potential build- 
up inside the mesoscopic conductor. It is the considera- 
tion of the internal potential distribution which has made 
the NEGF theory gauge invariant. At various limiting 
cases our theory predicts results consistent with those 
of scattering matrix theory and response theory, but our 
theory allows a general treatment of the full nonlinear DC 
transport regime which are, perhaps, impossible for the 
other formalisms. The present theory is natural to allow 
the inclusion of exchange and correlation potential within 
the density functional formalism. Hence it is appropriate 
for transport in multi-probe conductors in the quantum 
coherent regime, and we have applied it to the analysis 
of resonant tunneling with one and two resonance levels. 
Our results clearly show that without self-consistent anal- 
ysis the predicted current would depend on the choice of 
potential zero, which is wrong. For the tunneling device 
our analysis also indicated the importance of including 
charge polarization effect. This effect can be considered 
using the phenomenological model involving capacitance 
coefficients, as done here; or it can be included through 
numerical solutions of a charging model as carried out in 
Ref. pl^. Finally we found that in general a larger cur- 
rent is obtained when exchange and correlation effects 
are included into the analysis. 

Many further applications of the present formalism can 
be made. An important further development is to aban- 
don the wideband limit. In this work we have used this 
limit in order to derive analytical formula, but one can go 
beyond this limit in numerical calculations. In the wide- 
band limit, the coupling matrix is independent of energy. 
A consequence, as we observe from the I-V curves, is that 
the current increases monotonically with bias (no peaks 



in the I-V curve itself). Hence the negative differential 
resistance (NDR) can not be observed. To overcome this 
limitation, Jauho et. aJi3 have introduced a lower energy 
cutoff to allow a finite occupied bandwidth of the con- 
tact. This modification allowed NDR to appear, but the 
current at large voltage did not agree with experimental 
results. Therefore to obtain NDR quantitative correctly 
we must go beyond the wideband limit. This can be 
done within our formalism usipg_the numerical method 
developed by McLennan et. alEju. 
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FIGURE CAPTIONS 

(1) The current versus the voltage for a resonant tun- 
neling structure with two resonant levels for self- 
consistent solution and non-self-consistent solution 
with two different voltage gauges. Solid line: non- 
self-consistent solution for Vi = V and V2 = 0; dot- 
ted line: non-self-consistent solution for Vi = V/2 
and V2 = —V/2; dot-dashed line: self-consistent so- 
lution. The system parameters are Fi = T2 =0.1, 
El = 0.3, E2 = 1.4, and Ep = 0.11. 

(2) The corresponding differential conductance for the 
non-self-consistent solutions of Fig.(l). Solid line: 
non-self-consistent solution for Vi = V and V2 = 0; 
dotted line: non-self-consistent solution for Vi — 

V/2 and V2 = -V/2. 

(3) The gauge invariant current versus the voltage for 
a resonant tunneling structure with one resonant 
level. Solid line {Vxc not included) and dotted line 
(T4c included): Fi = F2 = 0.5, Ci = C2 = 5.0, 
AE = -1.0, a = 0.7; dot-dashed line {V^c not 



included) and dashed line {Vxc included): Fi =0.1, 
Fa = 0.5, Ci = 1.0, C2 = 5.0, AE = -1.0, a = 0.7. 

Fig. (4) The differential conductance of a resonant tunnel- 
ing structure with one resonant level. The param- 
eters are the same as that of Fig. (3). 

Fig. (5) The gauge invariant current versus the voltage for 
a resonant tunneling structure with two resonant 
levels. Solid line (Vxc not included) and dotted line 
(14c included): Fi = F2 = 0.1, Ci = C2 = 1.0, 
El = 0.3, E2 = 1.4, Ef = 0.11. 

Fig. (6) The differential conductance of a resonant tunnel- 
ing structure with two resonant levels. The param- 
eters are the same as that of Fig. (5). Inset: the ex- 
change and correlation potential Vxc versus voltage 
of a resonant tunneling structure with two resonant 
levels. 

Fig. (7) The differential conductance of a resonant tun- 
neling structure with Hubbard U term. The pa- 
rameters are: Fi = r2 = 0.5, Ci = C2 = 5.0, 
AE = -1.0, and Ui = 4.0. 
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